#include "centralized_polygon_system.h"

extern Config_parameters::Opti_alg Config_parameters::opti_algo;
extern double Config_parameters::gravity_scale;

namespace Field_packing
{
	const int Centralized_polygon_system::GRID_DIM = 100;
	const int Centralized_polygon_system::EXPONENT = 4;

	Centralized_polygon_system::Centralized_polygon_system(const std::vector<Packing_polygon*>& polygon_set, const Container_box& bndbox, double rm_scale, bool up_gravity)
		: System_potential_energy(bndbox),
		  dfg(bndbox, GRID_DIM, GRID_DIM),
		  //adf(polygon_set, 100),
		  nb_polygons(polygon_set.size()),
		  con_btm(bndbox.y_min()),
		  con_top(bndbox.y_max()),
		  areas(nb_polygons),
		  upward_gravity(up_gravity)
	{
		// grid field construction
		for (size_t i = 0; i < nb_polygons; i++)
		{
			dfg.append_distance_field(*polygon_set[i]);
		}

		// kd-tree field construction 
		
		// polygon sampling
		sample_pnts.resize(nb_polygons);
		centroids.reserve(nb_polygons);
		std::valarray<double> perimeters(0.0, nb_polygons);
		for (size_t i = 0; i < nb_polygons; i++)
		{
			for (size_t j = 0; j < polygon_set[i]->size(); j++)
			{
				perimeters[i] += CGAL::sqrt(polygon_set[i]->edge(j).squared_length());
			}
			areas[i] = polygon_set[i]->volume();
			centroids.push_back(polygon_set[i]->centroid());
		}
		const unsigned int sample_nb = 100; // how many points on the polygon with maximal area(approximately)
		rm = perimeters.sum() / nb_polygons / 100;
		rm *= rm_scale;
		//std::cout<<">>> r_m = "<<rm<<std::endl;
		//rm /= 8*sample_nb;
		double max_area = areas.max();
		double bd_square_length = std::max(bndbox.x_max() - bndbox.x_min(), bndbox.y_max() - bndbox.y_min());
		for (size_t i = 0; i < nb_polygons; i++)
		{
			Polygon_sampler ps(*polygon_set[i]);
			//ps.inner_sample(areas[i]/max_area*sample_nb, sample_pnts[i]);
			ps.inner_sample(bd_square_length/100, sample_pnts[i]);
		}

		// set up the optimization context
		int nb_constraints = 0;
		for (size_t i = 0; i < nb_polygons; i++)
			nb_constraints += polygon_set[i]->size();
		nb_constraints *= 2; // constraint on the x and y coordinate of vertices
		ktr_context = boost::shared_ptr<Knitro_context>(new Knitro_context(nb_polygons*3, nb_constraints, nb_constraints*3, 0));
		ktr_context->objType = KTR_OBJTYPE_GENERAL;
		ktr_context->objGoal = KTR_OBJGOAL_MINIMIZE;
		for (size_t i = 0; i < nb_polygons; i++)
		{
			ktr_context->xLoBnds[3*i] = -KTR_INFBOUND;
			ktr_context->xUpBnds[3*i] = KTR_INFBOUND;
			ktr_context->xLoBnds[3*i+1] = ktr_context->xLoBnds[3*i+2] = -KTR_INFBOUND;
			ktr_context->xUpBnds[3*i+1] = ktr_context->xUpBnds[3*i+2] = KTR_INFBOUND;
		}
		int *ptrCType = ktr_context->cType;
		double *ptrCLoBnds = ktr_context->cLoBnds, *ptrCUpBnds = ktr_context->cUpBnds;
		for (size_t i = 0; i < nb_polygons; i++)
		{
			for (size_t j = 0; j < polygon_set[i]->size(); j++)
			{
				*ptrCType++ = KTR_CONTYPE_GENERAL;
				*ptrCLoBnds++ = bounding_box.x_min();
				*ptrCUpBnds++ = bounding_box.x_max();
				*ptrCType++ = KTR_CONTYPE_GENERAL;
				*ptrCLoBnds++ = bounding_box.y_min();
				*ptrCUpBnds++ = bounding_box.y_max();
			}
		}
		int constraint_no = 0;
		int *ptrCons = ktr_context->jacIndexCons, *ptrVars = ktr_context->jacIndexVars;
		for (size_t i = 0; i < nb_polygons; i++)
		{
			for (size_t j = 0; j < polygon_set[i]->size(); j++)
			{
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i;
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i+1;
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i+2;
				constraint_no++;
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i;
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i+1;
				*ptrCons++ = constraint_no;
				*ptrVars++ = 3*i+2;
				constraint_no++;
			}
		}
		// polygon vertex for constraints
		//std::vector< Vector_2 > init_trans;
		//init_trans.reserve(nb_polygons);
		polygon_vtx.resize(nb_polygons);
		prev_trans.reserve(nb_polygons);
		// centralize the polygons
		// set the starting points for initialization
		for (size_t i = 0; i < nb_polygons; i++)
		{
			Point_2 c = polygon_set[i]->centroid();
			polygon_vtx[i].reserve(polygon_set[i]->size());
			for (size_t j = 0; j < polygon_set[i]->size(); j++)
			{
				polygon_vtx[i].push_back( polygon_set[i]->vertex(j) + Vector_2(c, CGAL::ORIGIN) );
			}
			ktr_context->xInitials[i*3] = 0.0;
			ktr_context->xInitials[i*3 + 1] = c.x();
			ktr_context->xInitials[i*3 + 2] = c.y();
			prev_trans.push_back(matrix<double>(2,2));
			matrix<double>& m = prev_trans.back();
			m(0,0) = 1.0;
			m(1,0) = 0.0;
			m(0,1) = c.x();
			m(1,1) = c.y();
		}
#ifdef _BOUNDARY_ENERGY_
		// sampling on the side boundary
		Point_2 lb(bndbox.x_min(), bndbox.y_min()), lt(bndbox.x_min(), bndbox.y_max()), 
				rb(bndbox.x_max(), bndbox.y_min()), rt(bndbox.x_max(), bndbox.y_max());
		size_t N = 20;
		bnd_pnts.reserve(2*(N+1));
		for (size_t i = 0; i <= N; i++)
		{
			double w = static_cast<double>(i)/N;
			Point_2 p = CGAL::ORIGIN + w*Vector_2(CGAL::ORIGIN, lb) + (1-w)*Vector_2(CGAL::ORIGIN, lt);
			bnd_pnts.push_back(p);
			p = CGAL::ORIGIN + w*Vector_2(CGAL::ORIGIN, rb) + (1-w)*Vector_2(CGAL::ORIGIN, rt);
			bnd_pnts.push_back(p);
		}
#endif
		// if this is just for simulated annealing, ignore gravity
		if (rm_scale > 2.0)
		{
			f = g;
			g = 0.0;
		}
	}

	void Centralized_polygon_system::move_to_new_point(const double *const x)
	{
		for (size_t i = 0; i < nb_polygons; i++)
		{
			double s = std::sin(x[3*i]), c = std::cos(x[3*i]);
			matrix<double>& m = prev_trans[i];
			Transformation_2 to_orig(CGAL::TRANSLATION, Vector_2(-m(0,1), -m(1,1)));
			Transformation_2 back_rot_orig(CGAL::ROTATION, -m(1,0), m(0,0));
			Transformation_2 rot(CGAL::ROTATION, s, c);
			Transformation_2 forward(CGAL::TRANSLATION, Vector_2(x[3*i+1], x[3*i+2]));
			Transformation_2 t = forward*rot*back_rot_orig*to_orig;

			// update the field
			dfg.update_distance_field(i, t);
			//adf.update_distance_field(i, t);

			// update the sampling points
			for (size_t j = 0; j < sample_pnts[i].size(); j++)
				sample_pnts[i][j] = t(sample_pnts[i][j]);
			// record the most recent transformation
			m(0,0) = c;
			m(1,0) = s;
			m(0,1) = x[3*i+1];
			m(1,1) = x[3*i+2];
			centroids[i] = t(centroids[i]);
		}	
	}
	double Centralized_polygon_system::molecular_energy_at(const double *const x)
	{
		//move_to_new_point(x);
#ifdef _CILK_
		cilk::reducer_opadd<double>  LJ_energy(0.0);
		cilk_for (size_t i = 0; i < nb_polygons; i++)
#else
		double LJ_energy = 0.0;
		for (size_t i = 0; i < nb_polygons; i++) // for every polygon
#endif
		{
			for (size_t j = 0; j < nb_polygons; j++) // for every sample point on each polygon
			{
				if (i == j) continue;
				for (size_t k = 0; k < sample_pnts[i].size(); k++) // compute potential energy with the other polygons pairwise
				{
					// grid based
					double d = dfg.distance_query(sample_pnts[i][k], j);
					if (d <= 0)	d = rm/50;
					
					// kd-tree based
					//Point_2 np;
					//double d = adf.nearest_point_at(sample_pnts[i][k], j, np);
					//d = std::sqrt(d);
					
					double rm6_over_d6 = std::pow(rm/d, EXPONENT/*6*/);
					LJ_energy += rm6_over_d6*(rm6_over_d6 - 2);
				}
			}
			//LJ_energy += areas[i]*(centroids[i].y() - con_btm)*g;
		}
#ifdef _CILK_
		return LJ_energy.get_value()*epsilon;
#else
		return LJ_energy*epsilon;
#endif
	}
	double Centralized_polygon_system::gravity_energy_at(const double *const x)
	{
		double gravity_energy = 0.0;
		for (size_t i = 0; i < nb_polygons; i++)
		{
			gravity_energy += areas[i]*(centroids[i].y() - con_btm);
		}
		return gravity_energy*g*Config_parameters::gravity_scale;
	}
	double Centralized_polygon_system::upward_gravity_energy_at(const double *const x)
	{
		double gravity_energy = 0.0;
		for (size_t i = 0; i < nb_polygons; i++)
		{
			gravity_energy += areas[i]*(con_top - centroids[i].y());
		}
		return gravity_energy*g*Config_parameters::gravity_scale;		
	}
	double Centralized_polygon_system::boundary_energy_at(const double *const x)
	{
#ifdef _CILK_
		cilk::reducer_opadd<double> bnd_energy(0.0);
		cilk_for (size_t i = 0; i < bnd_pnts.size(); i++)
#else
		double bnd_energy = 0.0;
		for (size_t i = 0; i < bnd_pnts.size(); i++)
#endif
		{
			for (size_t j = 0; j < nb_polygons; j++)
			{
				double d = dfg.distance_query(bnd_pnts[i], j);
				//if (d <= rm*2)	d = rm;
				if (Geex::Numeric::is_nan(d)) system("pause");
				double rm6_over_d6 = std::pow(rm/d, EXPONENT/*6*/);
				//std::cout<< rm6_over_d6 << std::endl;
				bnd_energy += rm6_over_d6*(rm6_over_d6 - 2);
			}
		}
#ifdef _CILK_
		//std::cout<<"boundary energy = " << bnd_energy.get_value() << std::endl;
		return bnd_energy.get_value()*epsilon;
#else
		return bnd_energy*epsilon;
#endif
	}
	void Centralized_polygon_system::boundary_energy_gradient_at(const double *const x, double *gradient)
	{
//#ifdef _CILK_
//		cilk_for (size_t i = 0; i < nb_polygons; i++)
//#else
		for (size_t i = 0; i < bnd_pnts.size(); i++)
//#endif
		{
			for (size_t j = 0; j < nb_polygons; j++)
			{
				double sj = std::sin(x[3*j]), cj = std::cos(x[3*j]);
				const matrix<double>& mj = prev_trans[j];
				Transformation_2 backj = Transformation_2(CGAL::ROTATION, -mj(1,0), mj(0,0))*Transformation_2(CGAL::TRANSLATION, Vector_2(-mj(0,1), -mj(1,1)));				
				Vector_2 grad;
				Point_2 closest_pnt;

				// grid based
				double d = dfg.gradient_query(bnd_pnts[i], j, grad, closest_pnt);
				//if (d <= rm*2)	d = rm;

				//double d = adf.nearest_point_at(sample_pnts[i][k], j, closest_pnt);
				//d = std::sqrt(d);

				double common_factor;

				double rm6_over_d6 = std::pow(rm/d, EXPONENT/*6*/);
				common_factor = -/*6*/EXPONENT*rm6_over_d6*(rm6_over_d6 - 1);
				common_factor /= d*d;
				common_factor *= epsilon;

				grad = Vector_2(closest_pnt, bnd_pnts[i]);
				grad = 2*grad;

				const Point_2& p = bnd_pnts[i];
				closest_pnt = backj(closest_pnt);
				gradient[3*j] += grad*Vector_2(closest_pnt.x()*sj+closest_pnt.y()*cj, -closest_pnt.x()*cj+closest_pnt.y()*sj)*common_factor;
				gradient[3*j+1] += -grad.x()*common_factor;
				gradient[3*j+2] += -grad.y()*common_factor;
			}
		}
	}
	void Centralized_polygon_system::molecular_energy_gradient_at(const double *const x, double *gradient)
	{
		//move_to_new_point(x);
#ifdef _CILK_
		boost::mutex mt;
		cilk_for (size_t i = 0; i < nb_polygons; i++)
#else
		for (size_t i = 0; i < nb_polygons; i++)
#endif
		{
			const matrix<double>& mi = prev_trans[i];
			Transformation_2 backi = Transformation_2(CGAL::ROTATION, -mi(1,0), mi(0,0))*Transformation_2(CGAL::TRANSLATION, Vector_2(-mi(0,1), -mi(1,1)));
			std::vector<Point_2> original_pnts(sample_pnts[i].size());
			std::transform(sample_pnts[i].begin(), sample_pnts[i].end(), original_pnts.begin(), backi);
			double si = std::sin(x[3*i]), ci = std::cos(x[3*i]);
			for (size_t j = 0; j < nb_polygons; j++)
			{
				if (i==j) continue;
				double sj = std::sin(x[3*j]), cj = std::cos(x[3*j]);
				const matrix<double>& mj = prev_trans[j];
				Transformation_2 backj = Transformation_2(CGAL::ROTATION, -mj(1,0), mj(0,0))*Transformation_2(CGAL::TRANSLATION, Vector_2(-mj(0,1), -mj(1,1)));				
				for (size_t k = 0; k < sample_pnts[i].size(); k++)
				{
					Vector_2 grad;
					Point_2 closest_pnt;

					// grid based
					double d = dfg.gradient_query(sample_pnts[i][k], j, grad, closest_pnt);
					if (d <= 0)	d = rm/50;

					//double d = adf.nearest_point_at(sample_pnts[i][k], j, closest_pnt);
					//d = std::sqrt(d);

					double common_factor;

					double rm6_over_d6 = std::pow(rm/d, EXPONENT/*6*/);
					common_factor = -/*6*/EXPONENT*rm6_over_d6*(rm6_over_d6 - 1);
					common_factor /= d*d;
					common_factor *= epsilon;

					grad = Vector_2(closest_pnt, sample_pnts[i][k]);
					grad = 2*grad;
					//if (d <= 0)
					//	g = -g;
					const Point_2& p = original_pnts[k];
					closest_pnt = backj(closest_pnt);
					boost::lock_guard<boost::mutex> lkg(mt);
					gradient[3*i] += grad*Vector_2(-p.x()*si-p.y()*ci, p.x()*ci-p.y()*si)*common_factor;
					gradient[3*i+1] += grad.x()*common_factor;
					gradient[3*i+2] += grad.y()*common_factor;
					gradient[3*j] += -grad*Vector_2(-closest_pnt.x()*sj-closest_pnt.y()*cj, closest_pnt.x()*cj-closest_pnt.y()*sj)*common_factor;
					gradient[3*j+1] += -grad.x()*common_factor;
					gradient[3*j+2] += -grad.y()*common_factor;
				}
			}
		}
	}
	void Centralized_polygon_system::gravity_energy_gradient_at(const double *const x, double *gradient)
	{
#ifdef _CILK_
		cilk_for (size_t i = 0; i < nb_polygons; i++)
#else
		for (size_t i = 0; i < nb_polygons; i++)
#endif
		{
			const matrix<double>& mi = prev_trans[i];
			Transformation_2 backi = Transformation_2(CGAL::ROTATION, -mi(1,0), mi(0,0))*Transformation_2(CGAL::TRANSLATION, Vector_2(-mi(0,1), -mi(1,1)));
			//std::vector<Point_2> original_pnts(sample_pnts[i].size());
			//std::transform(sample_pnts[i].begin(), sample_pnts[i].end(), original_pnts.begin(), backi);
			double si = std::sin(x[3*i]), ci = std::cos(x[3*i]);
			Point_2 orig_centroid = backi(centroids[i]);
			gradient[3*i] += areas[i]*g*(orig_centroid.x()*ci - orig_centroid.y()*si)*Config_parameters::gravity_scale;
			gradient[3*i+2] += areas[i]*g*Config_parameters::gravity_scale;
		}
	}
	void Centralized_polygon_system::energy_gradient_at(const double *const x, double *gradient)
	{
		std::fill(gradient, gradient + ktr_context->n, 0.0);
		move_to_new_point(x);
		molecular_energy_gradient_at(x, gradient);
		//if (g != 0.0)
		//{
			//std::cout<<"With gravity\n";
		if (!upward_gravity)
			gravity_energy_gradient_at(x, gradient);
		else
			upward_gravity_energy_gradient_at(x, gradient);
		//}
#ifdef _BOUNDARY_ENERGY_
		if (g != 0.0)
			boundary_energy_gradient_at(x, gradient);
#endif
	}
	void Centralized_polygon_system::upward_gravity_energy_gradient_at(const double *const x, double *gradient)
	{
#ifdef _CILK_
		cilk_for (size_t i = 0; i < nb_polygons; i++)
#else
		for (size_t i = 0; i < nb_polygons; i++)
#endif
		{
			const matrix<double>& mi = prev_trans[i];
			Transformation_2 backi = Transformation_2(CGAL::ROTATION, -mi(1,0), mi(0,0))*Transformation_2(CGAL::TRANSLATION, Vector_2(-mi(0,1), -mi(1,1)));
			//std::vector<Point_2> original_pnts(sample_pnts[i].size());
			//std::transform(sample_pnts[i].begin(), sample_pnts[i].end(), original_pnts.begin(), backi);
			double si = std::sin(x[3*i]), ci = std::cos(x[3*i]);
			Point_2 orig_centroid = backi(centroids[i]);
			gradient[3*i] += areas[i]*g*(-orig_centroid.x()*ci + orig_centroid.y()*si)*Config_parameters::gravity_scale;
			gradient[3*i+2] -= areas[i]*g*Config_parameters::gravity_scale;
		}	
	}
	void Centralized_polygon_system::constraint_at(const double *const x, double *constraint_values)
	{
		for (size_t i = 0; i < nb_polygons; i++)
		{
			double s = std::sin(x[3*i]), c = std::cos(x[3*i]);
			Transformation_2 t = Transformation_2(CGAL::TRANSLATION, Vector_2(x[3*i+1], x[3*i+2]))*Transformation_2(CGAL::ROTATION, s, c);
			for (size_t j = 0; j < polygon_vtx[i].size(); j++)
			{
				Point_2 p = t(polygon_vtx[i][j]);
				*constraint_values++ = p.x();
				*constraint_values++ = p.y();
			}
		}
	}
	void Centralized_polygon_system::constraint_grad_at(const double *const x, double *constraint_grad)
	{
		for (size_t i = 0; i < nb_polygons; i++)
		{
			double s = std::sin(x[3*i]), c = std::cos(x[3*i]);
			//Transformation_2 t(c, -s, x[3*i+1], s, c, x[3*i+2]);
			for (size_t j = 0; j < polygon_vtx[i].size(); j++)
			{
				const Point_2& p = polygon_vtx[i][j];
				*constraint_grad++ = -p.x()*s - p.y()*c;
				*constraint_grad++ = 1.0;
				*constraint_grad++ = 0.0;
				*constraint_grad++ = p.x()*c-p.y()*s;
				*constraint_grad++ = 0.0;
 				*constraint_grad++ = 1.0;
			}
		}
	}
	int Centralized_polygon_system::optimize(double *trans, bool apply_constraints)
	{
		KTR_context* kc = ktr_context->kc;
		if (kc == NULL)
			return -1000 ; // probably a license issue
		extern Config_parameters::Opti_alg opti_alg;
		if (KTR_set_int_param_by_name(kc, "algorithm", Config_parameters::opti_algo) != 0)
			return -1000;
		/* set options: automatic gradient and hessian matrix */
		if (KTR_set_int_param_by_name(kc, "bar_murule", 3) != 0) // "3" is plausible up to now if the above algo para is "1". "6" a bit less plausible
			return -1000;
		if (KTR_set_int_param_by_name(kc, "gradopt", KTR_GRADOPT_EXACT) != 0)
			return -1000 ;
		if (KTR_set_int_param_by_name(kc, "hessopt", KTR_HESSOPT_BFGS) != 0)
			return -1000 ;
		//if (KTR_set_int_param_by_name(kc, "bar_directinterval", 0) != 0)
		//	return -1000;
		if (KTR_set_int_param_by_name(kc, "bar_feasible", 1) != 0)
			return -1000;
		if (KTR_set_int_param_by_name(kc, "maxcgit", 5) != 0)
			return -1000;
		if (KTR_set_int_param_by_name(kc, "xtol", 1e-3) != 0)
			return -1000;
		// debug parameters
		if (KTR_set_int_param_by_name(kc, "debug", 0) != 0)
			return -1000;
		if (KTR_set_int_param_by_name(kc, "outlev", 0) != 0)
			return -1000 ;
		//if (KTR_set_int_param_by_name(kc, "outmode", 0) != 0)
		//	return -1000;
		//if (KTR_set_int_param_by_name(kc, "newpoint", 3) != 0)
		//	return -1000;
		if (KTR_set_int_param_by_name(kc, "bar_maxbacktrack", 10) != 0)
			return -1000;
		//if (KTR_set_int_param_by_name(kc, "bar_directinterval", 5) != 0)
		//	return -1000;
		int max_iter = 1000;
		if (g == 0.0) 
			max_iter = 3000;
		if (KTR_set_int_param_by_name(kc, "maxit", max_iter) != 0)
			return -1000 ;
		if (KTR_set_double_param_by_name(kc, "opttol", 1e-2) != 0)
			return -1000;
		//if (KTR_set_int_param_by_name(kc, "ms_enable", 1) != 0)
		//	return -1000;
		//if (KTR_set_int_param_by_name(kc, "ms_maxsolves", 5) != 0)
		//	return -1000;
		//if (KTR_set_int_param_by_name(kc, "ms_maxbndrange", 1) != 0)
		//	return -1000;
		/* register the callback function */
		if (KTR_set_func_callback(kc, &callback) != 0)
			return -1000 ;
		if (KTR_set_grad_callback(kc, &callback) != 0)
			return -1000 ;
		//if (KTR_set_newpoint_callback(kc, &newpoint_callback) != 0)
		//	return -1000;

		int nStatus = KTR_init_problem(	kc, ktr_context->n, ktr_context->objGoal, ktr_context->objType, 
										ktr_context->xLoBnds, ktr_context->xUpBnds, ktr_context->m,
										ktr_context->cType, ktr_context->cLoBnds, ktr_context->cUpBnds,
										ktr_context->nnzJ, ktr_context->jacIndexVars, ktr_context->jacIndexCons,
										ktr_context->nnzH, ktr_context->hessIndexRows, ktr_context->hessIndexCols, ktr_context->xInitials, NULL);
		if (nStatus)	return nStatus;
		nStatus = KTR_solve(kc, trans, ktr_context->lambda, 0, &ktr_context->obj, NULL, NULL, NULL, NULL, NULL, this);
		/* direvative check */
		//double init_obj = system_energy_at(trans);
		//std::fill(trans, trans + ktr_context->n, 0.0);
		//nStatus = KTR_check_first_ders(kc, trans, 1, 1.0e-10, 1.0e-10, 0, init_obj, NULL, NULL, NULL, this);
		/*************/
		// adjust gravity acceleration
		// current gravity potential
		if (g == 1.0)
		{
			double gravity_potential = 0.0;
			for (size_t i = 0; i < nb_polygons; i++)
			{
				gravity_potential += areas[i]*(centroids[i].y() - con_btm);
			}
			g = (ktr_context->obj - gravity_potential*g)/gravity_potential;
			g = std::fabs(g);
		}
		else if (g == 0.0)
			g = f;
		std::cout<<"current weight of gravity potential = "<<g<<std::endl;
		return nStatus;
	}

	void Centralized_polygon_system::draw_sample_pnts() const
	{
		glColor3f(1.0f, 0.0f, 0.0f);
		glPointSize(2.0f);
		glBegin(GL_POINTS);
		for (size_t i = 0; i < sample_pnts.size(); i++)
		{
			for (size_t j = 0; j < sample_pnts[i].size(); j++)
			{
				glVertex2d(sample_pnts[i][j].x(), sample_pnts[i][j].y());
			}
		}
		glEnd();
	}
}